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Abstract 

Using a generalized Rubinstein-Duke model we prove rigorously that kinematic disorder leaves 
the prediction of standard reptation theory for the scaling of the diffusion constant in the limit 
for long polymer chains D oc L~ 2 unaffected. Based on an analytical calculation as well as Monte 
Carlo simulations we predict kinematic disorder to affect the center of mass diffusion constant of 
an entangled polymer in the limit for long chains by the same factor as single particle diffusion in 
a random barrier model. 
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The derivation of large scale properties of polymer systems such as viscosity and diffusion 
constant out of microscopic models is among the basic problems of polymer science. For 
entangled polymers in gels or melts, de Gennesi predicted scaling laws for the dependence on 
the polymer length L for viscosity (rj oc L 3 ) and diffusion constant (D oc L~ 2 ). These scaling 
laws are assumed to be valid for the limiting case of polymer length going to infinity. In 
experiments, the apparent scaling laws D oc L~ 2A and 77 oc L 3A are foun d 2 ' 3 ' 4 . However, the 
experimental findings do not contradict the predictions of reptation, as a crossover due to 
decreasing finite size effects with increasing chain length cannot be ruled out 3 -. Presumably, 
contour length fluctuations (CLF) are one cause for the deviating scaling exponent for finite 
chain length 5 -^. In the framework of the repton model introduced by Rubinstein^ and 
Duk o 10 ' 11 (further on to be called RD model), which incorporates CLF, it is possible to 
calculate viscosity, diffusion coefficient and other quantities of interest. Good agreement both 
with theoretical and experimental results is foun d 6 ' 7 ' 14 ' 15 . Using the RD model with periodic 
boundary conditions Kooiman and van Leeuwen^ analytically calculated the proportionality 
constant c for the diffusion constant in the limit for infinite chain length: lim^oo DL 2 = c. 
They found c = w/(2d + 1), d being the dimensionality of the entanglement network and w 
a model constant defining the time scale. Building upon this result, Prahofer and Spohni 3 - 
rigorously derived the leading order term in 1/L for the diffusion constant in the RD model 
and furthermore proposed a scaling for the finite size effects: DL 2 — W/(2d + 1) oc L^ 13 , 
where 1/2 < (3 < 1. All results mentioned above were obtained under the assumption that 
the entanglement network which is topologically restricting the polymer is regular and static. 
Real entanglement networks such as gels or melts are disordered. The effects of a disordered 
environment can be manifold: 

• Spatial variations of the mobility of the 'defects' of stored length. 

• Locally fluctuating potential energy due to interactions between chain and environ- 
ment. 

• Entropically favorable regions of low entanglement density. 

• Relaxation of the environment (constraint release). 

Numerical investigations 16 showed that entropically favorable regions can for short chains 
substantially lower the diffusion constant by the creation of 'entropic traps'. So far, conclu- 
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sive investigations of diffusing polymers long enough to span several such traps are missing. 
Constraint release is considered to be of minor importance in gels but needs to be self 
consistently taken into account in melts 17 . Schafer, Baumgartner and Ebert^ numerically 
investigated the effect of kinematic disorder, i.e., disorder which affects the mobility of the 
chain segments while leaving the equilibrium distribution of chain configurations unaltered. 
Their investigation shows that reptation prevails in presence of kinematic disorder. How- 
ever, due to being based on MC simulations and thus relatively short chains no quantitative 
results could be obtained about the modifications the disorder would cause to the reptation 
prediction for the diffusion constant in the long chain limit. As reptation is shown to prevail, 
the scaling lim^oo-DL 2 = c must remain valid. The aim of this article is to calculate the 
coefficient c for a polymer diffusing in a disordered environment exhibiting kinematic disor- 
der. This means that we are interested in the behavior of the chain in the long chain length 
limit originally envisaged by standard reptation theory. Some of the results presented here 
in detail were briefly reported in an earlier communication^. 

The article is divided into six chapters: In chapter I the definition of the model is given and 
the master equation for the chain dynamics is presented in terms of the quantum Hamiltonian 
formalism's. This master equation yields the stationary distribution of chain conformations 
as shown in chapter II. The model with periodic boundary conditions is analyzed in chapter 
III. Adapting the strategy in Ref. [l2 for obtaining the diffusion constant for the periodic 
system yields a lower bound on the diffusion constant for the original system. Chapter IV is 
devoted to deriving an upp er bound to the diffusion constant by generalizing a variational 
approach used in Ref. Il3j for the ordered system. In chapter V Monte Carlo simulations 
of the polymer diffusion in kinematically disordered environments are presented. The last 
chapter addresses the dynamics of internal chain segments. 

I. A LATTICE MODEL FOR REPTATION WITH KINEMATIC DISORDER 

The RD model is a discretized model for reptation, i.e., it describes the dynamics of a 
polymer in an entanglement network. This network is assumed to be regular and static. 
In 3 dimensions it has the shape of a cubic lattice, where the polymer is forbidden to cut 
through the edges of the cubes. The faces of the cubes can by penetrated by the chain (see 
figure Q). The polymer itself is assumed to consist of L + 1 'reptons', i.e., segments of about 
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the entanglement length, which equals the lattice constant of the cubic lattice. 
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Figure 1: Two dimensional representation of a kinematically disordered lattice with a reptating 
polymer chain and the lattice gas mapping. Arrows indicate possible moves of the reptons. B = 
exp(E/2). 

The model defines dynamical rules for the reptons' motion. The dynamics on a smaller 
scale, i.e. high frequency Rouse modes, is averaged out. In this sense, the RD model is an 
effective model. The dynamical rules for the reptons are: 

a) Each cell occupied by the chain must contain at least one repton to ensure connectivity 
of the chain, since the length represented by a repton equals the lattice constant of the 
cells. The sequence of occupied pores corresponds to the tube as in standard reptation 
theory. 

b) End reptons can move to adjacent cells provided rule a) is not violated. 

c) Interior reptons can move to cells occupied by the neighboring reptons if allowed by 
rule a). This ensures the dynamics to be reptation, as any movements in the interior 
requires multiple reptons to occupy the same cell, which means that there is some 
excess 'stored length' available. This corresponds to the motion of 'defects' in the 
original de Gennes reptation model. 
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All moves are assumed to be thermally activated. In our model the cubic lattice repre- 
senting the entanglement network gets kinematically disordered by assigning an individual 
energy barrier to every boundary between adjacent cells. These barriers have to be overcome 
by a crossing repton. Note that we do not assign different energies to the reptons when re- 
siding in the cells. This choice of quenched disorder guarantees that only the mobility of the 
reptons, but not the equilibrium configurations are affected, as will be seen below. Let there 
be a E N different possible hopping rates W a , which are distributed randomly and occur 
with a probability f(W a ) throughout the lattice. We demand that the disorder averages 
(l/W) = Yla=i f{W a )/W a and (1/W 2 ) (the mean hopping time and its second moment) 
are finite. 

The RD model is used to describe the dynamics of an entangled polymer chain under 
the influence of an external electric field. A common example is DNA under electrophoresis 
conditions where the reptons carry a charge each and develop a drift velocity along the 
direction of an applied electric field. Let the field be oriented along the (111) diagonal 
of the cubes. We denote the (dimensionless) energy gain of a repton when moving from 
one cell to another along the direction of the field as E. By local detailed balance, moves 
across a cell boundary with assigned hopping rate W a along the field happen with rate 
W a exp(J5/2), those in the opposite direction with rate W a exp(— E/2)^. By projecting the 
reptons' positions on the direction of the field, their relative coordinates can be denoted as 
a 1-d lattice gas with L sites by the following prescription: 

• If the projected link between adjacent reptons i and i + 1 is oriented along (against) 
the field direction across a cell boundary with assigned rate W a , site i of the lattice 
gas is assigned the value yi = a (yi = —a). We interpret values ±a as particles. 

• If adjacent reptons i and i + 1 occupy the same cell, i.e., their projected positions 
coincide, the link is represented in the lattice gas as = at site i. We interpret a 
site, which is assigned a as being unoccupied. 

Thus the L + 1 coordinates of the reptons in direction of the field are translated into 
the equivalent set of the relative coordinates manifest as the assignment of a, —a or 
to the L sites of the lattice gas plus the center of mass coordinate's component in field 
direction. The dynamics of the reptons translates into the lattice gas picture as follows: In 
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the bulk particles of sort '±a' hop to the left with rate W a exp(±E/2) and to the right 
with rate W a exp(=fE/2), where each site can be occupied by at most one particle. The 
end dynamics in the lattice gas picture needs some care: Assuming y\ (y L ) to be non 
zero, the only possible move is the retraction of the end repton to the cell occupied by 
it's neighbor (rule a)). This retraction, being a particle annihilation event in the lattice 
gas picture, happens with the same rate as the respective move in the bulk. Assuming yi 
(ul) to be zero, the end repton can, according to rule b), move to any of the 2d adjacent 
cells. For half of these the move leads to links being along the field direction, the other half 
against it. The probability of the chosen move leading to the crossing of a cell boundary 
with rate W a being assigned to it is f(W a ). Thus the move of the repton, being a particle 
creation event in the lattice gas picture leads to y\ (y£) changing from '0' to '±a' with 
rate exp(=F E / '2) / (W a )W a d (exp(±E/2)f(W a )W a d). This choice of boundary dynamics 
is on average correct, but neglects the actual local structure of the network. In order to 
verify that the model is correct on long time scales (as relevant for the CMS diffusion) 
we performed Monte Carlo simulations comparing the diffusion constant of our projected 
model with d = 2 at different disorder distributions with a polymer chain moving according 
to the dynamics of the repton model in a two dimensional lattice with random but fixed 
hopping rates assigned to each cell boundary. It turns out that although the results differ 
for small chains (L < 10) they coincide within the statistical errors for longer chains. This 
legitimates our choice of boundary dynamics for investigating the behavior of long chains. 
Figure El shows an example where a disorder distribution with < W >= 0.505 was chosen. 

Every move of a particle in the lattice gas leads to a change in the component of the 
center of mass coordinate of the repton chain along the field direction. 

• Particle type 'a' moving to the right (left) decreases (increases) x by 1/(L + 1), as this 
is equivalent to a repton moving downward (upward). As there are L + l reptons each 
contributes 1/(L + 1) to the center of mass position. 

• Particle type '—a' moving to the right (left) increases (decreases) x by 1/(1/ + 1). 

In the subsequent chapters we will calculate the drift velocity v of the center of mass 
coordinate in presence of an applied electric field and by employing the Einstein relation the 
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Figure 2: Comparison between the diffusion coefficients for simulations of the disordered RD and 
two dimensional lattice model. 

zero field diffusion constant. When calculating v we restrict ourselves to the linear response 
regime, ignoring higher order field dependences. In order to calculate v only the change 
in the center of mass coordinate along the field x has to be known, but not the absolute 
value itself. Whereas the latter cannot be known from the lattice gas, the former is given 
by the difference of the currents of positive particles j + and that of the negative ones j~: 
v = j~ — j + . As the choice of the field direction is arbitrary, by use of the Einstein relation 
the model allows for computing the zero field diffusion constant along a distinct direction 
in rf-dimensional space. As zero field diffusion is isotropic this immediately yields the d- 
dimensional diffusion constant. This is due to distinguishing between particles of types a 
and —a in contrast to the original projected repton model used by Rubinstein, which allows 
only for the computation of the curvilinear diffusion constant along the contour of the chain. 
In the Rubinstein model additional assumptions are necessary to relate the curvilinear to the 
rf-dimensional diffusion constant. Our model as well as the RD model allows for computation 
of the latter quantity within the model. 

II. QUANTUM HAMILTONIAN AND STATIONARY STATE 

The model introduced in the previous chapter describes a Markov process and thus the 
dynamics can be written in form of a master equation. For convenience we will use the 
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quantum Hamiltonian formalism 20 to write down the master equation at zero electric field 
and solve for the stationary state. 

At each site of the lattice gas with length L, 2a + 1 values for are possible. Therefore, 
the state space X has the dimension (2a + 1) L . Every rj £ X is assigned a vector \i] > 
and a transposed vector < 77 1 . These vectors constitute a basis for the space of system 
configurations X. A probability distribution P(rj) = P v can thus be written as a probability 
vector: 

\P>=J2 P ( r l)\v> ■ (1) 

Let a summation vector < s\ be defined as < s\ = ^2 veX < ^1- Normalization of the 
probability vector is given if < s\P >= 1. The generator for the dynamics of the system 
is the matrix H: The off-diagonal elements of the matrix contain the (negative) transition 
rates t(rj,rj') from a state rj to rj: 

<r } \H\rf>=H Virf = -t(ri,rf). (2) 
The diagonal elements H n ^ contain the sum of all outgoing rates : 

<r ] \H\r l >=H ri:11 = J2 t (v',v)- (3) 

It is easily checked that conservation of probability, i.e. < s\H = in the language of 
the quantum Hamiltonian formalism, is fulfilled by H. The master equation 

J t P v(t)= E Uv,v')PAt)-Kv\v)Pv(t)) (4) 
describing the Markovian dynamics of the system can be written as: 

j t \P(t) >= -H\P(t) > . (5) 
The stationary state \P* > is thus characterized by the equation 

H\P* >= 0. (6) 
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In this formalism the expectation value < F(t) >= ^ E(i])P(i], t) °f an operator F 
is written as follows 20 : The operator F : X — > X is represented by a diagonal matrix 
F = ^ n F{rj)\r] >< t)\. Then 

< F(t) > =< s\Fe~ m \P(0) > 

(7) 

=< s\e m Fe- Ht \P(0) > . 
Let the time-dependent operator F(t) for t > be defined as 

F(t) = e Ht Fe- Ht . (8) 

< > is an expectation value that it is not only averaged over possible realizations of 
the process but also over the initial states according to P{0). In the following it is assumed 
that |P(0) >= \P* >, so that 

< F(t) >=< s\F(t)\P* > . (9) 

For our model, we choose a tensor product basis as follows: Let at each site of the lattice 
gas the unit vector e± denote yi = 0, t2a denote yi = a and e2 a +i denote = —a. A 
state vector for a state rj = (1, — 3, .., 5, — 2) for example can then be written as \r) >= 
z<i ® e-j ® ... ® eio <8> e$). Using this basis, the following operator creates a particle of type a 
at site i, provided it was previously unoccupied: 

oi iX (l) = l®^S(2a,l) = E {2aA) (i), (10) 

i—l L—i 

where E(2 a ,i) is the matrix with a single entry 1 at row 2a and column 1. Similarly, 
the operator a„_i(i) = E(2 a +i,i)(i) creates a particle of type —a, if possible. The corre- 
sponding annihilation operators at site i are a a> i(i) = En^^i) and a Q) _x(z) = E^ 12a +i)(i)- 
To formulate the diagonal part of the quantum Hamiltonian the matrices u(i) = E^^(i), 
v a ,i{i) = E(2a,2a)(i) an d v a -i(i) = -E(2a+i,2a+i) (i) are employed. Thus, the quantum Hamil- 
tonian for the model defined in the previous chapter at zero field reads: 
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a 
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+W a (-a a>s (L) +v ai8 (L))] 
Plugging a product measure ansatz 
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into -Popenl-Popen(O) >= leads to a very simple set of equations for the probabilities p a ,±i 
and one finds 



P* (0) >= 

open\ I 



( 1 \ 

df{W x 



XL 



{2d+l) L ' 



(13) 



df(W a ) 
\df(W a )J 

The geometrical equilibrium conformation of the chain depends on the probability of 
occurrence for links between reptons along or against the field, irrespective of the assigned 
hopping rates of possibly crossed cell boundaries. This means we have to consider the overall 
probability for particles of positive sign at a site, which is Yla=i f(W<x)d/ (2d+l) = dj (2d+l), 
or negative sign, being also d/(2d+ 1), respectively. These are equal to the probabilities 
found for the disorder free RD-model 22 which shows that the chosen kind of disorder is indeed 
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kinematic disorder, as it influences only the mobility but not the equilibrium configurations 
of the chain. 

Note that it is not possible to compute the stationary state for the model with a non-zero 
field at arbitrary L. Therefore the Einstein relation 1/L d/dE(< v(E) >) = D is not a 
feasible way of straightforwardly computing the diffusion constant. Extra input is needed 
to overcome this problem. 

III. MODEL WITH PERIODIC BOUNDARY CONDITIONS 

In this section we adapt a calculation due to Kooiman and van Leeuwen^ for the original 
RD-model with periodic boundary conditions to the case with kinematic disorder. The 
quantum Hamiltonian for the periodic system at zero field reads 

a L 

H per = ^2^2 W a[-ai )S (i + l)a a ^('i) - al s {i)a a>s (i + 1) 

s=±l a=l i=l 

+u{i + l)v a>s (i) + u(i)v a , s {i + 1)] . (14) 

Plugging |-P o * pe „(0) > as for the open system into the equation H per \P* pen (0) >= for 
the periodic system shows that at zero electric field the stationary state of the open system 
is also stationary with respect to the dynamics of a corresponding periodic system. In a 
periodic system, the phase space is non ergodic, as neither the order nor the number of 
occurring particles on such a ring can be changed. Therefore, every connected subset of 
the phase space ('channel') has its own stationary state. We can calculate the stationary 
state for the periodic system in presence of a field E by mapping the system to a disordered 
zero range process, as introduced by Benjamini, Ferrari and Landim 23 . This means that 
the Einstein relation can be employed to obtain the diffusion constant in the periodic case. 
Here, the definition of the drift velocity of the center of mass and the corresponding diffusion 
constant are induced from the system with open boundary conditions: v = j~ — j + . 
Instead of characterizing the system by the spins on the lattice gas y = (yi, ■■,yt) it can 
equivalently be characterized by the sets s = (si, .., sl) containing the signs of the non zero 
Hi, w = (wi, ..,wl) containing the absolute values of these yi and n = (ni, ..,n^) where rii 
amounts to the number of '0' between Sj and Sj + i on the lattice gas. Thus every lattice 
gas configuration on a ring can be translated into a configuration on a lattice of length M 
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in the zero range picture: Site % carries rii zero range particles and is separated from site 
i + 1 by a bond characterized by Sj+i and Wi + \. The total number K of zero range particles 



must equal the total number L — M of '0' in the lattice gas: YldLi ni = K = L — M. The 
dynamics of the lattice gas picture translates into the zero range picture as follows: The 
configuration (.., rij, rij+i, ..) changes to (..,rij — 1,^+1 + 1, ••) with rate hJ^Wj+i and to 
(..,rij + 1,^+1 — 1, ..) with rate hj+iWj+i, where hj = exp(—Esj/2). This means that the 
random hopping rates as well as the Sj are not assigned to individual particles, but to bonds 
between sites in the zero range (zr) picture. At E — the zr-particles move as in a random 
barrier energy landscape. Moves of the zr particles cause changes of the center of mass 
coordinate as defined above: A zr-particle hopping to the right across a bond with Sj > 
(sj < 0) increases (decreases) the center of mass position by 1/L. Conversely, a zr-particle 
hopping to the left across a bond with Sj > (sj < 0) decreases (increases) the center of 
mass position by 1/L. As with the field free case of the open system, we use a quantum 
Hamiltonian and a tensor basis for the state space to compute the stationary state of the zr 
lattice of length M with a given s and w. Let e±(i) denote an unoccupied site % and e m (i) 
where m > 1 a site i occupied by rii = m — 1 particles. Here, vectors are infinite-dimensional. 
The matrix for creation of a particle at site i is then given by: 



b] = 



... 

1 ... 



\ 







1 



(15) 



The matrix for annihilation of a particle at site i is accordingly given by: 



A) 1 ... ^ 



bi = 



1 







(16) 
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For constructing the diagonal part of the Hamiltonian we need the following type of 
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matrix (note that < s\b\ =< s\l and < s\bi =< s|toj) 

fn n 



rrii 



0.. 

o i • 



••. ... 
1 '•• 



(17) 



A hopping event from site i to site i + 1 is described by the combined action of the 
matrices hi and b\ +1 This yields the following expression for the quantum Hamiltonian: 



= ^ ( - hi +1 b i+1 b\w i+1 - h^hbl^Wi+i + h^Wi+Wi + h i+1 w i+ im i+1 j 

i=l ^ ' 



A yet unnormalized product ansatz 



yields 



IP* >= 

-L. 



M 














® . . . <g> 










w 







(18) 



(19) 



M 



H \P* > 

11 zr\ 1 zr ^ 



i=l 



mi|P z * r > -hJiW. 



Zi+l 



m i+1 \P* r > 



+h i l l w i+l m i \Pl r > +hi +1 w i+1 m i+1 \P* r > ^ 

M 

= E 



i=l 



rrii\P zr > -h i Wi rrii\P zr > 



+/i i+ 1 1 u> i+ im i |P* r > +h, i w, i m i \P* r > ^. 



\P* r > is stationary provided 

h i+1 w i+1 z i+1 + h^WiZi-i = h^Wi+iZi + hiWiZi 



(20) 



(21) 
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Applying a general solution 24 to the case at hand and including normalization the sta- 
tionary state for a channel characterized by s and w and with a total number of particles 
K yields 



M 



\P* zr ( S ,w,K,M)>= f[z? i \n 1 ,...,n M >— 1 - , (22) 

(ni,...,n M ) i=1 ^(n njf) Hi=l Z i 



where 



M i-1 

z, = Y^]^T—Jl h 'Uj- (23) 

i=l n l+i w ^ j=1 

The primed sums are meant to be summations under the constraint Yli=i rii = K. Know- 
ing the stationary state for each channel, the drift velocity for the individual channels can 
be computed. In the lattice gas picture, every time a particle Sj = —1 (s$ = 1) is hopping 
to the right (left) it changes the position of the center of mass by 1/(L + 1). The opposite 
process, i.e. a particle Sj = —1 (sj = 1) hopping to the left (right) changes the center of 
mass position by — 1/(L+ 1). Therefore, as mentioned above, v = j~ i.e. the difference 
between the currents of negative and positive particles. This operator v translates to the zr 
picture as follows: 



J ~ 3 + -> 7 Yl s *+i w *+i( /l *+i 6 A t +i ~ h i+1 b i+1 bl). (24) 



4 = 1 

Thus, the drift velocity of a given channel is given by 



< v(s,w, K,M) > = < s\-Y,( s i+i h i+i w i+i b ib\ +1 

-s i+1 h i+1 w i+1 b i+1 b])e(K)\P: r > — (25) 

Ck 

Here, Q(K) projects on those states of \P* r > which have a constant number of particles 
K and ck is the normalization. So Q(K) / ck\Pz T >= \Pzr( s i w > K) >• Using this form allows 
us to make use of the fact that due to the combined effect of b;b\ which is to redistribute 
the particles without changing their number, we can commute b^- and Q(K). This allows 
us then to apply the matrices to the product measure, where its effect is easy to see. 
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< v(s, w,K,M) > = < s\^-^2(s i+1 h i ^ 1 w i+1 b i bl +1 

i=i 

-s i+ ihi + iWi +1 bi + ibl)Q(K)\P* r > — 
1 M 

= < s \j^2( s i+i h i+iWi+i®( K )bib\ +1 



L 

i=l 

-Si + ih i+1 Wi +1 <d(K)b i+1 b\)\P* r > — 

ck 

I 

T22( S i+l h i+l w i+l^- < S \ m i+l\®i K ) P zr > 
Li . z i+l 

-s i+1 h i+1 w i+1 ^ < s\ mi \Q{K)P: r >)— (26) 



We now have to calculate the expression ck for the product state: 

M 

c *= E IR ( 2? ) 

(ni,...,n M ) i=l 

where the primed sum again means summing with the constraint YlfLi n i = K- 
A similar type of summation is found when regarding the sums of the type 

< S |m^|P; r >= Y. ^IK- («) 



(m,...,n M ) 3 =1 



where the double primed sum has the constraints YljLi n i = K an d n /3 7^ 0- A straight- 
forward calculation of the sums is impossible due to the constraints. To simplify the task, 
we can profitably use the following identity^: 



m , M 



(ni,...,ra M ) 1=1 (ni,...,n M ) j=l 

Thus we can transform the sums into integrals: 
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(ni,...,n M ),W/3^0 J=l 

1 1 f da J^r 1 

7=1 



j=l 

Qk-\,m hn 

= • (31) 

The integrals Qk,m satisfy the recursion relation 

Qk,m = Qk,m-i + zmQk-i,m- (32) 
An explicit solution of this relation yields: 

M M , 

Qk,m = Y. z ? +M ~ 1 n^-^) -1 - (33) 
1=1 1=1 

The primed product indicates i ^ I. 
From this point on we are not going to carry through the complete calculation for 

< v(s, w, K, M) >. We will expand it into a series in E and keep only first order terms 
which is sufficient for employing the Einstein relation. We start from the expression for 

< v(s, w, K, M) > after having inserted the result for 



I ts\ 1 i -1 Zi Qk-i,m 

<v{s,w,K)> = - > (s i+1 h i+1 w i+1 z i+l — 



Zi+i Qk-1,M 
—Si + illi + iWi + i Zi 



Zi Qk,m 
1 \" (u-l u ^Qk- I W 

= j ) J s i+1 (h i+1 w i+1 Zi - h i+1 w i+1 z i+1 ) 



T / j Vi+1 "'i-h±~J ■"i-t-±>~i-t-±~i-t-±,/ ^ 



1=1 

M 



M 



Q 



K-l.M 



1 

= 7 J> + i(l-exp(-£S)) 
^ t=i 

= f(l-ex P (-^))%^£, (34) 

where S = J^fLi s i- F° r the second equality we used the explicit form of the z t . The 
term in brackets is easily expanded into a series in E: 1 — exp(ES') = ES + o(E 2 ). We 
keep only the first order term in E. This means that when we are expanding the expression 
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Qk-i,m/Qk,m we need only to keep zero order terms in E, as all other contributions will 
vanish when using the Einstein relation. The only terms containing E in the Qk,m are the 
Zi, for which we find to first order in E: 



m j-l 

,2 



Y J 



j=l h i+j w i+i fc= l 
M . „ j-l 



j=i Wl+ i k=i 

M 

+o(E). 

/ — ' w i+j 



3=1 

= z + o{E). (35) 

Here, z is defined as the zero field value of Z{. The fact that all Zi are equal at zero field 
is decisive for the explicit integration of Qk,m- We obtain: 

M 



3=1 



M + K-l 
M-l 



z K . 



(36) 



It is now clear how to proceed with evaluating < v(s, w, K, M) >, as 

n (K+M-2\ M-l „ 1 

QK-1,M _ V M-l J Z _ 1 

This yields for < v(s, w, K, M) >: 



(37) 



<v(s,w,K,M)> = \ Y,s 1+1 ES k+ K m _ i \ + o{E^) 

i=l 

- ^nbr;^)- (38) 

This is our final expression for < w(s, w,K, M) >. Knowing the drift velocity for each 
channel characterized by s,w we have to give each of those a weight \l/(s, w, K, M) when 
averaging over the channels. Following an argument by Prahofer— , we show below that 
provided a configuration in the zr picture is weighted such that at E = its probability as 
contained in P r * ng (s, w, K, M) = ^(s, w, K, M)P* r (s, w, K, M) is equal to the probability 
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of the corresponding state in the lattice gas of the open system as given by Po pen (0), for the 
diffusion constants the following relation holds: -D^+i ^ D p ^ r . 

When relating zr to lattice gas probabilities it has to be kept in mind that the periodic zr 
system has M sites, while the corresponding periodic lattice gas has L > M sites. This 
means a state of the zr system with probability q zr can be permutated M times where 
due to translational invariance all resulting states have the same probability. A lattice gas 
configuration with probability qi g can be cycled through L permutations leading to equal 
probabilities. Therefore, equal probabilities of respective configurations means: Mq zr = 
Lqig. The weight factor thus has to be chosen as follows: 

*( S .w,/^fl=p^( M )n/W. (39) 
Thus we find for the drift velocity when averaged over the channels: 



v = J2 Yl *(s,w) <v(s,w, K,M) > 

M s=(si,...,s M ) w=(w 1 ,...,w M ) 

v- 1 ^., i M (L\L-M /Vtt „ a\ 

M v / \ / \ w j— i J 



(40) 



This can be rewritten as 



L 1 

v = E ^(M)M(-). (41) 



m=i ' z 



Here Q(M) contains all the factors depending on M occurring in the previous expression 
and < 1/ z > is a disorder average. f2(M) is a function that is sharply peaked at M = L/2 
implying that in the limit for large L only terms with large M significantly contribute to the 
result. In the limiting case of L — > oo invoking the central limit theorem (remember that 
we demanded < 1/W > and < 1/W 2 > to be finite) yields: 

(z) = M(i/wy (42) 

Therefore 



A 1 ^ r d M ( L\L — M 1 1 
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Using the Einstein relation yields in the limit for large L 



11 11 
m = L(L-l) (2rf+l) (1 "(2rf+l)^ ) <J) 

= (44) 
\w) 

where D* is the diffusion constant for the ordered cas o 12 4 3 . 
Thus we have shown that Iuxyl-,^ D L 2 = l/[(2d+ 1) < 1/W >} for the periodic case, which 
is a lower bound to the open case. Note that the factor 1/ < 1/W > is the same as occurring 
in the single particle diffusion constant for the random barrier model. In this model random 
energy barriers are assigned to bonds between sites just as in the zr picture of our problem. 

IV. MODEL WITH OPEN BOUNDARY CONDITIONS 
A. Variational formula 

In order to find an upper bound on the diffusion constant we follow the strategy of Ref. 



13l . The quantum Hamiltonian H open for our model in its representation using the tensor 
basis can be decomposed into the sum of a diagonal part (D) and a non diagonal one (M). 
Each of these can themselves be split into a part describing those moves leading to an 
increase of the center of mass coordinate (M + and D + ) and a respective part connected to 
a decrease (M~ and D~). Thus H = D+ + D~ - M+ - M~ . 
The following statement? 2 ^ holds for the diffusion constant D: 



D = inf 



<D + + D- > - 2 < s\[D + - D~]g\P* > + < s\gHg\P* > 



_2(L+1) 2 ' ~ (L + l) 



inf F[g]. 

gen 



(45) 



Here, Q is the space of diagonal matrices with dimension (2a + 1) L . Thus plugging in any 
diagonal matrix g into the functional F[g] yields an upper bound on D. The challenge is to 
choose g such that the upper bound gets as small as possible. Due to the dimensionality of 
the diagonal matrices increasing exponentially with L, an exact minimization is not feasible. 
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Still, some information can be gained by observing that the functional is convex and thus 
the minimizing go G Q is the only matrix for which SF[g] vanishes. Choosing go such that 



< s\g H =< s\[D + - D-] (46) 

the variation vanishes. Unfortunately, this formula cannot be solved for go- We will show 
below how still some insight can be gained from this equation. Introducing the matrix d, 
where d(y'\y) =< y'\d\y > gives the change in the CMS coordinate when a transition from 
state y to y' is made and the matrix w, where w(y'\y) denotes the corresponding rate, the 
variational formula may be written as: 

D= inf i £ w(y'\y)P*(y)[d(y'\y) + g(y')-g(y)] 2 , (47) 

y',yex 

where g(y) =< s\g\y >. Following Ref. [22 we now prove that £)° pen > D p ^_ v To simplify 



the notation we remark that in our model with periodic boundary conditions as in the 
original RD model transitions from a given state y are only allowed to a state y' = y l,l+1 , 
where the spins yi and y i+ i are interchanged. Thus the last formula applied to the periodic 
case reads as 



dt +1 = mf \j2[Y, w ^ l+x \y) p ^y^ ( 48 ) 

1=0 \yeX J 

Here, each transition between states in a lattice gas of length L + 1 changes the cm 
coordinate by ±1/ (L+ 1), as there are as many reptons as bonds between them. In contrast 
to that in the open boundary case there is one repton more as there are bonds leading to the 
cm coordinate changing by ±1/ (L + 1) in any transition of a lattice gas of length L. Here, 
transitions between states can not only occur by interchanging spins y\ and y%+\ in the bulk, 
but also by creation and annihilation events at the ends. These events can equivalently be 
viewed as interchanges of the spins y\ or y L respectively with 'imaginary' spins y and yL+i 
provided the rates for these interchanges at the ends are adapted such that they correspond 
to the creation and annihilation rates as demanded by the model. Taking the distribution 
at the 'imaginary' sites as for any other site in the lattice gas (remember that we have a 
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homogeneous product measure) the introduction of a parameter 

{2d + 1 if 1 = or I — L 
(49) 
1 for Z = 1,..,L-1 

leads to a fulfillment of this requirement and the functional for the diffusion constant of 

the open chain of length L reads as 



D° L pen = mf ^ q ( w (y^ +1 |y)P*(y)[^(z/' i+1 |y) + 9(y l,l+1 ) ~ g(y)\ 2 j (so) 

Z=0 Vj/GX / 

where u> and d are as for the periodic case. Letting y = (y Q , .., i/l+i) for open and periodic 
system and observing that as y does not occur in the functional for D v ^ x the averaging of 
yo over P* yields 1 and thus does not change the result. 

D°r[g] - DTIM 
= co [J2 w (y 0,1 \y) p *(y)[d(y 0A \y) + g(y°> 1 ) -g(y)f 

\yex 

+c L (^ W (/' L+1 b)P*(y)[% L ' L+1 b) + ^(/' L+1 )-^(y)] s 

- £ ^(/ +1 ' 1 b)^(y)[% ' 1 b) + aiy ' 1 ) - g{y)f 
-J2 w (y L ' L+1 \y) p *(y)[d(y L > L+1 \y) +g(y L > L+1 ) -g(y)} 2 

y&X 

+^ £ w (y L,L+1 \y)P*(y)\d{y L)L+1 \y) + g(y L ' L+1 ) — g{ 

y&X 

> (51) 

Therefore, the diffusion constant of the periodic chain with length L + 1 is a lower bound 
for the open chain with length L. 

B. Functional for diffusion with kinematic disorder 

The space of diagonal matrices Q from which an appropriate g G Q has to be chosen 
such that F[g] gets as small as possible, has dimension (2a + 1) L . A scalar product on 
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can be denned for arbitrary m, n G Vt as < s\mn\P* >. A basis for Q can be built from the 
matrices yi(i) to y~2<r{i) acting non trivially only at site i and the unit matrix. We specify 
only the following: 

Va{i) = ^-^E (2a)2a) (i) - -l—E {7a+li7a+1) (i) for a = 1, ..,a, (52) 



2d „ 2 ^ 1 



2/(7+1 W 



2d+ 1 v ' n ' ^ 2d+ 1 

J =2 



For the remaining a — 1 matrices we demand that they are chosen such that for all 
a,b— 1, ..,2cr , i,j — 1, -,L and a ^ b we have < y a (i)yb{j) >= as we U as < 2/a(*) >— 0. 
It is easily checked that the defined matrices are mutually orthogonal. A basis for the space 
of diagonal matrices Q is then given by: 

'2(7+1 ^ 

II%( 7 ^C {!,••,£} Vj;/«n//, = . (54) 

j=i J 

The choice of a # G fulfilling < s|[-D + - D~] =< s\gH leads to 5F[g] = 0. In spite of 
this relation not being solvable for g information can be gained by observing that Q contains 
subsets ^i 1; ..,j CT which are invariant to H in the sense that if g G f^..^, then the diagonal 
matrix g' fulfilling < s\gH =< s\g' is in f^,..,^ too, where 



2(7+1 . 

^u,..,v = s P an j II vAh^h c {!> ••>} Vj;/ Q n/ /3 = Va^/3;|I,-| =ij Vj = l,..,(jL 
^ j=i J 

(55) 

By calculating < s|[-D + — -D - ] we find that the minimizing g is in the subspace 
©J=i ^u,..,j CT where i& = 5^ V /c = l,..,cr. Thus the most general ansatz for g based 
on this information is: 

a L+l 

9 = J2Y1 £ ?Wi,.,Vi^( i )^+i(Ui)-y2 ff +i(^+i) (56) 

a=l i=0 7 ct+1 ,..,7 2ct +i 

where the coefficients c are real numbers and the sets I a , .., ha+i are mutually disjoint. 
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C. Ansatz for g 



The general ansatz given above for g contains too many parameters for letting a mini- 
mization appear feasible. Instead we use a generalization of the ansatz used by Prahofer— 
for the ordered RD model: 



cr , L+l L+l , 

j-E E a ^) + E 1 

a=l ^ i=0 fc,fc'=0 



ctk>y a {k)i)o+i(k') 



(57) 



L+l , 

+ E e o,p,g^(°)^+l(p)^+l(^) 
o,p,g=0;p<g 

The primed sums indicate the summation variables to be mutually different. Note that 
this ansatz reflects a particular choice of parameters for the general ansatz. Inserting the 
trial function into F[g] yields the following functional: 



Dig] 



E{(^t/«E 

a=l u ; Z=0 



+ 



L + l 

2d 



Q 

2d 



{ W IMA W. 



(2d + 



1)2 E" ( C " fc ' C ?+l,k> 



(2d+l) v ~ M+1 
2d 



(2d+l) 



e Z,(Z+l,fc') e Z+l,(Z,/c')- 



2tf 



+Ep5Ti)/(«wE-(«4.- 

/3=1 v 7 A: 



^ 1 



4d 



P<<? 



4d 2 



(58) 



+ E ?2dTi^ /W) E«a, ? ) - <^) 2 

,3=1 V ; o,<j 

The free parameters are the coefficients af,c% k , and e° which are, generalizing results 
in Ref. Q chosen as 



H+l 



E 

z=o 



2d 



L + l (2d + l) 



-1,1+1 H+1.L 



(59) 



a? = 
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(2d+l)C , 
2d 1 



i,3 



'L+l-i,L+l-j 



(2rf+l) 2 C // o \ (L+l-p)(L+l-q) 
idn y \n> (L+l) 2 (L-o) 



'L+l-o,L+l-q,L+l—p 



for i > j 

for o < p < q 
for p < q < o. 



(60) 



(61) 



Here, g is a monotonically decreasing smooth function on the real numbers with g(x) = 1 
for x < which is decreasing exponentially fast for x — ► oo. This function has the property 
that 



E 

i=0 



0(n) 



(62) 



where r, s are integers, and n = (L + 1) . Given this choice, the following upper bound 



for D°l™ is found 



proper i ^ 
u L+\ — 



1 



< w >, 



(63) 



(2d + 1) (L + l) 2 
where < > is the disorder average. 

This completes the proof that asymptotically D oc 1/L 2 in the presence of kinematic 
disorder. 



V. MC RESULTS 

In the case of the ordered RD model the results for D in the periodic case^ coincide to 
leading order with the upper bound for the open case^. In our model we managed to give 
bounds on D which both scale with 1/L 2 in the limit for long chains. This implies that 
kinematic disorder does not ruin the scaling relation of standard reptation theory. Still, 
the numerical prefactors differ. We believe that it is the result for the lower bound which 
correctly describes the diffusive behavior for long chains. This lower bound was obtained by 
a rigorous calculation while the upper bound results from a variational treatment. Given 
that our conjecture is true and for long chains D = l/[L 2 (2d + 1) < 1/W >} then for any 
choice of disorder distribution DL 2 < 1/W > plotted against the chain length is constant. 
We performed MC simulations with different disorder distributions. It turns out that L = 30 
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is sufficient to observe the curves all converging to the same constant value. For the cases 
< l/W >= 1, < l/W >= 0.02 and < l/W >= 0.03 figure El shows our MC results. These 
results suggest that in the long chain length limit D = l/[L 2 {2d + 1) < l/W >}. 



m 

j 

10 20 30 40 

L 

Figure 3: D(0)L 2 < l/W > versus L for the cases as described in the text. 



VI. DYNAMICS OF INTERNAL SEGMENTS 

The surprising result of the previous sections is the observation that the effect of kinematic 
disorder on the collective behavior of all connected polymer segments is (to leading order in 
system size) the same as on a simple point like object in the same disorder environment. In 
order to understand this observation we now consider the dynamics of the internal segments 
in the hydrodynamic limit of vanishing lattice spacing. For the local concentration of reptons 
one obtains from the usual ordered repton model in this limit Rouse dynamics^ restricted 
to motion inside the tube 5 . The boundary dynamics, i.e., the hopping into and out of the 
tube resp. describe the entropic tensile force acting on the chain ends and keeping the 
polymer in its stretched equilibrium conformation. For understanding the dependendence of 
the diffusion coefficient on kinematic disorder, which affects mostly the bulk of the polymer 
chain, it is sufficient to focus on the hydrodynamic behavior of the bulk reptons. In order 
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to study this limit it is convenient to first investigate the associated zero-range process and 
then translate the result into the hydrodynamic limit of the exclusion process. 

For the zero-range process described in Section III with E = 0, the average number of 
particles pi(t) at site i obeys the exact time evolution equation given by 

= W t z^(t) + W i+1 z i+1 (t) - (Wi + W i+1 ) Zl {t) (64) 

where Zj(t) is the probability that site i is occupied at time t. It is known that in the steady 
state, the occupancy probability zf s is spatially uniform and can be related to the steady 
state particle density as^ 

ss = _J\ — (65) 

1 i + pr 

It follows that the steady state density profile is uniform in spite of disorder. 

To understand the dynamics in the hydrodynamic limit of vanishing lattice spacing, we 
expand Eq. (j64|) to second order in lattice constant and find 



dp(x 1 t) d 
dt dx 



(66) 



dx 

At large enough time, the system is expected to be in local equilibrium so that we 
may assume the usual approximation for z(x,t) — p(x,t) in the steady state to be valid: 
z(x, t) « p(x,t)/(l + p(x,t)), where equation (|63|) was used. We are interested in the den- 
sity fluctuations about the steady state, Ap(x, t) = p(x, t) — p ss (x) where p ss (x) = p = K/M. 
Retaining the lowest nonvanishing term in the expansion of the preceding equation in powers 
of Ap(x, t), we obtain 



dAp(x, t) 1 d 



■ dA P (x,t) 
ox 



(67) 



dt (1 + p) 2 dx 

The above equation describes a random walker in one dimension, diffusing in a random 
medium with bond-symmetric hopping rates W(x). It can be shown that at large time and 
length scales, the random walker can be described by a single, effective diffusion constant, 
V = l/[(l + p) 2 (1/W)), provided (1/W) is finite 25 . Thus we obtain the bulk diffusion 
constant to be given by T> in the zero range process. 

Regarding the site as particle and mass as hole clusters, the above model maps onto 
symmetric exclusion process (SEP) with particlewise disorder. We want to calculate the 
bulk diffusion constant in the SEP picture using the above results for zero range process. 
Since the steady state density profile is uniform in both the pictures, the average local 
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density rij in the vicinity of the location of particle i in SEP is related to that in the zero 
range process as = 1/(1 + p,,). Then the density fluctuation about the steady state in the 
vicinity of particle % is given by 

An, = ^-^Ap 4 + 0((Ap) 2 ) . (68) 

We note that the above is true only at large enough times, as was pointed out in a similar 
analysis for the tagged particle correlation function for SEP without disorder—. 
Using Eq.flEZJ and Ea. ljBSjl . we find that 



dAn(x,t) 1 d 



■ dAn^ 
ox 



(69) 



dt (1 + p) 2 dx 

Further note that x is the space index in the zero range process, while it labels the particles 
in the SEP. The space coordinate y in symmetric exclusion process is related to x as 



2/ ~ / dx'p(x') + x , (70) 

which gives 

| = (1 + ^ + 0(AB) . (71) 

Thus the particle density fluctuations in the SEP obey 

dAn(y,t) d 



dt dy 



W(y) dAnM 



(72) 



dy 

Using the random walker analogy, we obtain the effective diffusion constant T> at large times 
to be equal to 1/ (1/W) in the symmetric exclusion process. 

This calculation shows that the internal segments of the polymer chain perform Rouse- 
type dynamics also in the presence of kinematic disorder, but with a disorder-dependent 
diffusion coefficient. This explains the occurrence of the same correction to the diffusion 
coefficient for the long time behavior of the polymer chain as a whole. 



VII. CONCLUSIONS 



It is the aim of this paper to disentangle the effects of the various types of disorder 
which one may expect to have significant impact on the dynamics of systems of entangled 
flexible polymers. We have focused on kinematic disorder which leaves the equilibrium 
conformation unchanged compared to a hypothetical ordered entanglement network (which 
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could in principle be manufactured artificially by placing a single polymer in an ordered 
array of obstacles on a surface). We have proved (rigorously in terms of the RD- model for 
reptation) that the asymptotic length scaling of the diffusion coefficient of the polymer chain 
remains as predicted by standard reptation theory. By studying the hydrodynamic limit we 
have shown that the individual polymer segments inside the tube perform Rouse dynamics in 
a disordered environment which corresponds to a system of local random barriers. Therefore 
the amplitude of the diffusion coefficient becomes dependent on the disorder in the same 
way a single particle in a random barrier system. 
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